Critical behavior of self-assembled rigid rods on two-dimensional lattices: 
Bethe-Peierls approximation and Monte Carlo simulations 
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The critical behavior of adsorbed monomers that reversibly polymerize into linear chains with 
restricted orientations relative to the substrate has been studied. In the model considered here, 
which is known as self-assembled rigid rods (SAARs) model, the surface is represented by a two- 
dimensional lattice and a continuous orientational transition occurs as a function of temperature 
and coverage. The phase diagrams were obtained for the square, triangular and honeycomb lattices 
by means of Monte Carlo simulations and finite-size scaling analysis. The numerical results were 
compared with Bethe-Peierls analytical predictions about the orientational transition for the square 
and triangular lattices. The analysis of the phase diagrams, along with the behavior of the critical 
average rod lengths, showed that the critical properties of the model do not depend on the structure 
of the lattice at low temperatures (coverage), revealing a one-dimensional behavior in this regime. 
Finally, the universality class of the SAARs model, which has been subject of controversy, has been 
revisited. 
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I. INTRODUCTION 

Self-assembly has become a topic of increasing inter- 
est in recent years. One reason for this interest is that 
it is central to understanding structure formation in liv- 
ing systems Q. As a consequence, a significant research 
effort has been devoted to enhance our understanding 
of the theoretical basis of the fundamental mechanisms 
governing self-assembly and the observables required to 
characterize the interactions driving thermodynamic self- 
assembly transitions 0, Q . More related to the present 
work, several research groups reported on the assembly 
of particles in linear chains • Despite of the number 
of contributions to this problem, the knowledge of how 
this process works is still limited. 

It is obvious that a complete analysis of the self- 
assembly phenomenon is quite a difficult subject because 
of the complexity of the involved microscopic mecha- 
nisms. For this reason, the understanding of simple mod- 
els with increasing complexity might be a help and a 
guide to establish a general framework for the study of 
this kind of systems, and to stimulate the development 
of more sophisticated models which can be able to repro- 
duce concrete experimental situations. 
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In this context, an interesting model was introduced 
by Tavarcs ct al. @. The system in Ref. 0] consists of 
monomers with two attractive (sticky) poles that poly- 
merize reversibly into polydisperse chains and, at the 
same time, undergo a continuous isotropic-nematic (IN) 
phase transition. Using an approach in the spirit of the 
Zwanzig model , the authors studied the IN transition 
occurring in the system. The obtained results revealed 
that nematic ordering enhances bonding. In addition, 
the average rod length was described quantitatively in 
both phases, while the location of the ordering transi- 
tion, which was found to be continuous, was predicted 
semiquantitatively by the theory. Finally, Tavares et al. 
assumed as working hypothesis that the nature of the IN 
transition remains unchanged with respect to the case 
of monodisperse rigid rods on square lattices, where the 
transition is in the 2D Ising universality class j7l| - [l^ . 

From the seminal work by Tavares et al., a series of 
papers exploring the self-assembled rigid rods (SARRs) 
model have been published [l5l - [23j . These studies can 
be separated in two groups: (i) those dealing with the 
nature and universality of the phase transition occur- 
ring in the system |15l - [2l| and (ii) those dealing with 
the temperature-coverage phase diagram of the system 

up. 

With respect to the first point, the universality class of 
the model has been a subject of controversy. Thus, the 
criticality of the SARRs model in the square lattice was 
investigated in Ref. [l5| by means of canonical Monte 
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Carlo (MC) simulation and finite-size scaling (FSS) the- 
ory. The existence of a continuous phase transition was 
confirmed. In addition, the determination of the critical 
exponents, along with the behavior of the Binder cumu- 
lants, revealed that the universality class of the IN tran- 
sition changes from 2D Ising-type for monodisperse RRs 
without self-assembly to g = 1 Potts-type (random perco- 
lation) for polydisperse SARRs. Later, a multicanonical 
MC method based on a Wang-Landau sampling scheme 
was used by Almarza et al. [1^1 to reinvestigate the crit- 
ical behavior of the model studied in Refs. [1, Em- 
ploying the crossing point of the Binder cumulants and 
the value of the critical exponent of the correlation length 
(i/), it was observed that the criticality of the SARRs 
model is in the 2D Ising class, as in models of monodis- 
perse RRs [i3,[i|. 

The results in Refs. [HI, [3, along with the recent 
study in Ref. [l7j . indicate that the system under study 
represents an interesting case where the use of differ- 
ent statistical ensembles (canonical or grand canonical) 
leads to different and well-established universality classes 
{q = I Potts- type OT q = 2 Potts-type, respectively). A 
similar scheme was observed for triangular lattices, where 
canonical MC simulations indicated that the IN transi- 
tion of SARRs at intermediate density belongs to the 
q=l Potts universality class [l^ In contrast with 
this result, a g = 3 Potts- type universality was obtained 
by using grand canonical MC simulations plj . 

Among the studies of the second group, the 
temperature-coverage phase dia gram of SARRs on square 
lattices was calculated in Ref. By using MC sim- 

ulations, mean-field theory, and a renormalization group 
(RSRG) approach, the critical line which separates re- 
gions of isotropic and nematic stability was obtained and 
characterized. The results showed that the theory pre- 
sented in Ref. Q overestimates the critical temperature 
in all range of coverage. Small differences appear between 
simulation and theoretical results for small values of 6; 
however, the disagreement turns out to be significantly 
large for larger ^'s. On the other hand, RSRG reproduces 
qualitatively the shape of the critical line, but systemati- 
cally underestimates the critical temperature. The main 
prediction of RSRG approach is that the critical proper- 
tics of the whole line are associated to a unique second- 
order fixed point, confirming the continuous nature of the 
transition. Concerning the MF results, the theory pre- 
dicts the existence of a first-order transition line and a 
tricritical point. This finding is in sharp contrast to that 
obtained by MC simulations and RSRG approach. 

More recently, the main adsorption properties of 
SARRs on square and triangular lattices have been ad- 
dressed The study demonstrated that the adsorp- 
tion isotherms appear as sensitive quantities to the IN 
phase transition, allowing to reproduce the temperature- 
coverage phase diagram of the system for square lattices, 
and to obtain a first determination of the phase diagram 
for triangular lattices. 

Following the line of Refs. [l5l - [23| . the present pa- 



per deals with the two aspects above mentioned. On 
one hand, the problem of the universality is revisited, 
clarifying recent results obtained for triangular lattices 
(with conclusions that can be extrapolated to the hon- 
eycomb lattice case). On the other hand, the complete 
temperature-coverage phase diagrams were obtained for 
the square, triangular and honeycomb lattices by means 
of MC simulations and FSS analysis. The critical lines 
were also calculated for the square and triangular lat- 
tices within the Bethe-Peierls (BP) or quasi-chemical ap- 
proach, as formulated in the Cluster Variational Method. 
Comparisons with MF and RSRG data indicate that BP 
represents a qualitative advance in the analytical descrip- 
tion of the phase diagram of SARRs. 

The rest of the paper is organized as follows. In Sec. II 
we describe the lattice-gas model. The simulation scheme 
and computational results are given in Sec. III. In Sec. 
IV we present the analytical approximations, and com- 
pare the MC results with the theoretical calculations. In 
Sec. V we review previous results on the universality 
class of the SAARs model in the triangular and honey- 
comb lattices. Finally, the general conclusions are drawn 
in Sec. VI. 



II. LATTICE-GAS MODEL 

As in Refs. d, HB-dll, a system of self-assembled rods 
with a discrete number of orientations in two dimensions 
is considered. The substrate is represented by a square, 
triangular or honeycomb lattice oi M ~ L x L adsorp- 
tion sites, with periodic boundary conditions. N particles 
arc adsorbed on the substrate with m possible orienta- 
tions along the principal axes of the array, being to ~ 2 
for square lattices and to = 3 for triangular and honey- 
comb lattices. The rods interact with nearest-neighbors 
(NN) through anisotropic attractive interactions. Thus, 
a cluster or uninterrupted sequence of bonded particles 
is a self-assembled rod. Then, in the canonical ensemble 
the adsorbed phase is characterized by the Hamiltonian 

H = -wY^ {l% • • l}, (1) 

(id) 

where indicates a sum over NN sites; w represents 
the NN lateral interaction between two neighboring sites 
i and j; the energy is lowered by an amount \w\ only if 
the NN monomers are aligned with each other and with 
the intermolecular vector r^j; CTj is the occupation vec- 
tor, with (7j = if the site j is empty and cfj = Xk 
if the site j is occupied by a particle with orientation 
parallel to Xk- Xk are a set of unit vectors along the 
crystalline axes. In Eq.([T]) , div represents an integer di- 
vision, so the result inside { } can be either or 1 (i.e. 
the fractional part is discarded). The integer division is 
redundant in the case of the Hamiltonian for the square 
lattice, but it avoids additional lateral interactions [l^l 
that promote the condensation of the monomers in the 
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triangular and honeycomb lattice [20[ , restricting the at- 
tractive couplings only to those pairs of NN monomers 
whose orientations are aligned with each other and with 
the monomer-monomer lattice direction, in line with the 
model in the square lattice. 

The grand canonical Hamiltonian of the model is given 

by 



H = -w 



, • (7j | div l| - - Eo)^ |<Ti 



(2) 

where eo is the adsorption energy of an adparticlc on a 
site and ^ is the chemical potential. In the present work, 
€0 = and /i is the only parameter that determine the 
strength of the adsorption. 

On the other hand, it is worth mentioning that, while 
the concept of linear rod is trivial for square and trian- 
gular lattices, in a honeycomb lattice the geometry does 
not allow for the existence of a linear array of monomers. 
In this case, we call linear rod to a chain of adjacent 
monomers that can be assembled in only three types of 
sequences, defining three directions in similar way to the 
triangular lattice. For more details about the model in 
the honeycomb lattice, see Ref. [l9| . 

In the case of a square lattice the grand canonical 
Hamiltonian ^ can be exactly mapped into a spin-1 
model mi 



^ = -I E m + S,){S^^ + S,){x,.n,)+ 

<i,j> 



where Si = 0, ±1 and xi,X2 are unit vectors along the 
two orthogonal crystalline directions. 5^ ± 1 represent 
the vertical ((7^ = £2) and the horizontal {at = xi) ori- 
entations, while Si = represents the empty state. 

In the case of a triangular lattice the model can be 
formulated in terms of a diluted g = 3 anisotropic Potts 
model. We associate to each site of the lattice a spin 
variable <7i = 0,1,2,3, such that ai = represents the 
empty state and the states (Tj = 1,2,3 represent a bar 
oriented along the three unit vectors Xk [k ~ 1,2,3), 
where the angle between any pair of them is 27r/3 . Then 
the Hamiltonian can be written as 



H = -w y^<5(gi,g)(5(gj,o-)(5(rij-,±Xg) - 

<i.,j> (7=1 

-M^[l-<5(a„0)] (4) 



where S{a,a') is the Kronecker delta function. These 
alternative representations of the model are useful for 
the analytical treatment. 



III. MC SIMULATIONS 

A. MC method 

We have used a standard importance sampling MC 
method in the canonical and grand canonical ensembles 
[25} and FSS techniques [2^. All calculations were car- 
ried out using the parallel cluster BACO of Universidad 
Nacional de San Luis, Argentina. 

1. Canonical MC simulations 

MC simulations in the canonical ensemble were used 
to obtain the results presented in the subsection III.B. 
The procedure is as follows. Starting with a random ini- 
tial configuration (sites occupied with concentration 9 = 
N/M and particle axis orientation chosen at random), 
successive configurations are generated by attempting to 
move single particles (monomers) . One of the two (trans- 
lation or rotation) moves is chosen at random. In a 
translation move, an occupied site and an empty site are 
randomly selected and their coordinates are established. 
Then, an attempt is made to interchange its occup ancy 
state with probability given by the Metropolis rule j27l |: 
P = min {1, cxp (— /3AiJ)}, where AH is the difference 
between the Hamiltonians of the final and initial states 
and /3 = I/UbT (being ks the Boltzmann constant). For 
a rotation move, the rotational state of the chosen parti- 
cle is changed with a probability determined by Metropo- 
lis rule. A MC step (MCS) is achieved when 6 x M sites 
have been tested to change its occupancy state. Typi- 
cally, the equilibrium state can be well reproduced after 
discarding the first 5 x 10^ MCS. Then, the next 6 x 10^ 
MCS are used to compute averages. 

2. Grand canonical MC simulations 

MC simulations in the grand canonical ensemble have 
been carried out to understand the discrepancy between 
the results of [l^ and [2l| about the universality class 
of SAARs model in the triangular lattice (Sec. V). The 
procedure is as follows. For a given pair of values of 
T and ^, an initial configuration with N monomers ad- 
sorbed at random positions and orientations (on M sites) 
is generated. Then an adsorption-desorption process is 
started, where the lattice sites are tested to change its 
occupancy state with probability given by the Metropo- 



lis rule [27|: P = min {1, cxp (-/3Ai?)}, where Ai? is 



the difference between the Hamiltonians of the final and 
initial states and (3 = l/ksT. Insertion and removal of 
monomers, with a given orientation, are attempted with 
equal probability. For this purpose, an axis orientation 
(with probability 1/3 for the triangular lattice) and a 
lattice site are chosen at random. If the selected lattice 
site is empty, an attempt to place a monomer (with the 
orientation previously chosen) on the site is made. If, 
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instead, the site is occupied, then the algorithm checks 
the orientational state of the adsorbed monomer and if 
this coincides with the previously chosen orientation, an 
attempt to desorb the particle is performed; otherwise, 
the trial ends. A MCS is achieved when M sites have 
been tested to change its occupancy state. The equilib- 
rium state can be well reproduced after discarding the 
first 6 X 10^ MCS. Then, averages are taken over 6 x 10^ 
successive configurations. 



3. Order parameters, Binder cumulant and FSS analysis 

In order to follow the formation of the ncmatic phase 
from the isotropic phase, we use the order parameter de- 
fined in Ref. [l5| for the square lattice. 



\N.,-Nh\ 



(5) 



where Nh{Ny) is the number of monomers aligned along 
the horizontal (vertical) direction, and N is the number 
of total monomers on the lattice {N = Nh + N^). 

For the triangular and honeycomb lattices, the order 
parameter has been defined in Ref. jl^ as: 



Q 



\ni + n2 +n3\ 



n2 



l"3| 



(6) 



where each vector fii is associated to one of the 3 possi- 
ble orientations (or directions) for a chain on the lattice. 
In addition: 1) the n^'s lie in the same plane (or are co- 
planar) and point radially outward from a given point 
P which is defined as coordinate origin; 2) the angle be- 
tween two consecutive vectors. Hi and fti^i, is equal to 
27r/3; and 3) the magnitude of fii is equal to the num- 
ber of rods aligned along the z-direction (for a graphical 
representation, see Fig. 1(a) in Ref. [13 )• Note that 
the ni's have the same directions as the q vectors in Ref. 

In our canonical MC simulations, we fixed the density 
9, and monitored the order parameter (Q) as function 
of temperature T. The reduced fourth-order (Binder) 
cumulant Ul [2a|i "^^s calculated as: 



c/L(r) = i- 



3(Q2)2' 



(7) 



where the thermal average (...) means the usual time av- 
erage throughout the MC simulation. 

The critical behavior of the SAARs model has been 
investigated by means of FSS analysis. The FSS theory 
implies the following behavior for (Q) and Ul at critical- 
ity: 



and 



(Q)=L-^/-Q(Li/^e), 



(8) 



(9) 
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FIG. 1: Size dependence of the order parameter as a function 
of temperature for = 0.5 and 9 — 1 (inset). 



for L ^ oo, e — !• such that L^/^t~ finite, where e = 
T/Tc — 1 for canonical MC simulations and e = fi/ fic — ^ 
for grand canonical MC simulations. Here (3 and i/ are 
the standard critical exponents of the order parameter 
HQ) ^ (~£)'^ for £ ~^ 0~i L ~^ oo)i and correlation 
length ^ (^ ~ le]"" for e — > 0,L — !> oo), respectively. Q 
and Ul are scaling functions for the respective quantities. 

Finally, we calculated the average rod length on the 
transition line that, at fixed coverage, increases as the 
temperature decreases. At each MCS the average rod 
length may be calculated as 



Hmcs = 



N 



N - [Nb. 



>onds 



{L)rods\ 



(10) 



where N is the number of monomers adsorbed on the 
lattice; iV^onds is the number of bonds between pairs of 



nearest-neighbors monomers. N^ 



{L)rod& 



is the number of 



rods with length L; its inclusion prevents counting spu- 
rious bonds introduced by the periodic boundary con- 
ditions (i.e., bonds that do not contribute to the rod's 
length). Thus the equilibrium average rod length is ob- 
tained from 



i={eMcs), (11) 



B. Computational results 

1. Behavior of the Binder cumulant 

The critical behavior of the SAARs model has been 
investigated by means of the computational scheme de- 
scribed in the previous section for the canonical ensemble 
and FSS analysis [1^ [2^ . In order to illustrate the be- 
havior of the Binder cumulant in the critical regime, we 
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FIG. 2: Curves of Ul vs T* for = 0.5 (a) and 6* = 1 (b). 
From their intersections one obtained T* . In the insets, the 
data are plotted over a wider range of temperatures. 



of Ul{T*), which is independent of the system size at 
T* ~ T*. In other words, T* can be found from the in- 
tersection of the curve Ul{T*) for different values of L, 
since U* = Ul{T*) =const. As an example, Fig. 2 shows 
the reduced four-order cumulants U l plotted versus T* 
for the cases studied in Fig. 1. The values obtained 
for the critical temperature were T* = 0.220(2) (corre- 
sponding to 6* = 0.5) and T* = 0.476(3) (corresponding 
to 61 1). 

The behavior of the reduced fourth-order cumulant as 
a function of temperature also allows to make a prelim- 
inary identification of the order and universality class of 
the phase transition occurring in the system [25|. Thus, 
the curves in Fig. 2 exhibit the typical behavior of the 
cumulants in the presence of a continuous phase tran- 
sition. Namely, the order parameter cumulant shows a 
smooth drop from 2/3 to 0, instead of a characteristic 
deep ( neg ative) minimum, as in a first-order phase tran- 
sition j25j . 

The value of the intersection point U* shows two dif- 
ferent behaviors, which can be visualized from Fig. 2. 
On one hand, the value obtained for U* at — 0.5 
{U* = 0.643(4)) is consistent with the g = 1 Potts univer- 
sality class observed in Ref. [l^, where the system 
was studied at a fixed temperature (T* « 0.222). On the 
other hand, at = 1, the fixed value of the cumulants, 
U* = 0.605(6), is more consistent with previous estimates 
for the three-state Potts model (see for instance Ref. [lO] , 
where U* =0.613 [3l|). However, even though the value 
of U* may be taken as a first indication of universality, 
a detailed calculation of critical exponents is required 
for an accurate determination of the universality class. 
In Sec. V, the distinction between the two universality 
classes above is considered based on the determination of 
the critical exponent of the correlation length. 



show here the results for the triangular lattice case at 
intermediate and full coverage. In addition, these results 
will be useful when we address the question of whether 
or not the universality class can change as a result of 
the constant density constraint applied in the canonical 
ensemble (Sec. V). 

Figure 1 (inset of Fig. 1) shows the behavior of the 
order parameter versus the reduced temperature T* = 
kBT/w for several lattice sizes and 6 = 0.5 [l^ {6 = 1). 
As it can be observed, (Q) appears as a proper order 
parameter to elucidate the phase transition. When the 
system is disordered (T* > T* , being T* the critical tem- 
perature), all orientations are equivalent and (Q) is zero. 
In the critical regime (T* < T*), the particles align along 
one direction and {Q) is different from zero. 

Hereafter we discuss the behavior of the critical tem- 
perature as a function of coverage. The standard theory 
of FSS allows for various efficient routes to estimate T* 
from MC data [lH, [2^. One of these methods, which will 
be used in this case, is from the temperature dependence 



2. Phase diagrams of SAARs on square, triangular and 
honeycomb lattices 

The isotropic-nematic phase transition in a model of 
self-assembled rigid rods with restricted orientations was 
considered for the first time by Tavares et al. 01. The 
temperature-coverage phase diagram in Ref. |9|, ob- 
tained using an approach in the spirit of the Zwanzig 
model [loj . is qualitative only, and the theory overesti- 
mates the critical temperature in all ranges of coverage 
(especially at high coverages) [13 • However, for small 
values of 0, small differences appear between simulation 
and theoretical results (Fig. 3). As seen in Fig. 3, the 
critical lines separate regions of isotropic and nematic 
stability, and show that the nematic phase is stable at 
low temperatures and high densities. In addition, the 
phase diagrams show a marked increase of the density 
difference from dilute isotropic to dense nematic phases 
upon increasing the attraction between monomer units 
(i.e., decreasing the temperature). 

The critical line for SARRs on the square lattice, was 
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FIG. 3: Temperature-coverage diagrams for the SARRs model 
on different lattices. For comparative purposes, the critical 
curve reported by Tavares et al. 0] is shown as a solid line. 
The open and solid squares are from Refs. [S^] and [31, re- 
spectively. The open triangles and hexagons represent the 
results obtained in this work for the triangular and honey- 
comb lattices, respectively. Solid triangles are from Almarza 
et al. [2l|] (in black) and from Ref. [23 | (in grey). 



obtained by means of numerical simulations in Refs. 22 1 
and [l^. Fig. 3 shows the critical line reported in 22 1 
and only one point (the lowest coverage obtained) of the 
critical line reported in [l^. As explained in Ref. [ig} . 
the line of critical points continues beyond the lowest 
density showed in Fig. 3, however, the rapid increase of 
the average length of the rods at low densities and tem- 
peratures prevents an efficient simulation of the system. 

For the case of the triangular lattice, some critical 
points were obtained, at high coverages in Ref. [2l| and 
at intermediate coverages in Ref. (see Fig. 3). In the 
latter case, the points were obtained from the singular- 
ities in the adsorption isotherms. To corroborate these 
previous results and complete the phase diagram con- 
struction, the procedure used in III.B.l (to obtain the 
critical temperature) was repeated for 9 ranging between 
0.2 and 1. The same procedure was done for the case of 
the honeycomb lattice, this way the complete phase dia- 
gram of honeycomb lattice is reported here for the first 
time. All results are collected in Fig. 3. Together, the 
phase diagrams show that the critical properties coincide 
in the very low-temperature (coverage) regime. 



3. Critical average rod lengths 

Using a generalization of the theory of associating flu- 
ids, Tavares et al. @ obtained an analytic expression for 
the equilibrium average rod length in the SARRs model, 
with orientation a along the xi,X2 directions in two di- 
mensions (2D). In Fig. 4 (inset), we present a comparison 
between theoretical [9[ and numerical results for the equi- 
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Critical average rod lengths. See description in the 



fibrium average rod length on the transition line. In the 
numerical case, the points were obtained by using the 
Eqs. pU)) and ([TT|) and square lattices. By comparing 
these curves it can be seen a qualitative agreement up 
to 9 ~ 0.8. However, at higher coverages, MC simula- 
tions show a gradual increase of the critical average rod 
length, to a value of £ « 7 (at = 1). This value is 
interesting because it coincides with the minimum value 
of k {kmin = 7), which allows the formation of a nematic 
phase in long straight rigid rods of length k (fc-mers and 
monodisperse case) , on a square or triangular lattice [TTI - 

The average rod length on the critical line was also 
obtained for triangular and honeycomb lattices, see Fig. 
4. Two observations can be made from Fig. 4: (i) At 
intermediate coverage (9 w 0.5, dashed line), the crit- 
ical average rod lengths for the square and triangular 
lattices are similar and near to 7 (the minimum value 
of k which allows the formation of a nematic phase in 
monodisperse rigid rods, on a square or triangular lat- 
tice |lll - [l3| ). On the other hand, the critical average rod 
length (at 9 « 0.5) for the honeycomb lattice, is near to 
11, which coincides with the minimum value of k for the 
existence of a nematic phase in the case of monodisperse 
rigid rods on the honeycomb lattice |14| . (ii) Although 
the trend is more clear for the square and triangular lat- 
tices, all lines tend to converge as the coverage decreases 
towards zero; revealing a one-dimensional behavior in the 
very low-temperature (coverage) regime. As was shown 
in Ref. (2l| . at zero density limit, where the average rod 
length diverges, an equilibrium polymerization transition 
occurs. 
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IV. ANALYTICAL APPROXIMATIONS 

In this section we ealeulate the phase diagram within 
the Bethe-Peierls (BP) or quasichcmical approximation. 
To do that we use the Cluster Variation Method (CVM) 
[33| . In the CVM the BP approximation is obtained min- 
imizing a variational free energy expressed in terms of 
reduced probability densities, namely 

F = Trpff + fcBr|(l-gc)^Tr,pfhogpf' + 

+ E^-^i^°g^s|' (12) 

<i,j> J 

where Qc is the coordination number of the lattice, p[^^ 

(2) 

and pl j are one and two sites reduced densities respec- 
tively and it is assumed that p ~ Y[<i j> Pi^j ■ Pi^'' ^^.d 

(2) 

PI ■ can be expressed in terms of local one and two site 
averages, which are used as variational parameters and 
are related through the reducibility conditions: 



x^j = {SiSj) ^Tri,jSiSjpl^^, (17) 

y,, EE {SfS]) ^Tr,,, SfS^f^l (18) 

z,, EE {S,S])^Tn,,S.S^p^], (19) 

U, EE {S,Sf)=Tt,^,S,Sfp^^, (20) 

and imposing the reducibility conditions (|13p . the two 

(2) 

site density functions pl ■ can be obtained in terms of 
the variational parameters m, 6, Xij, yij, Zij and Uj 
(see Supplementary Material). Replacing into Eq. p2)) 
we obtain an expression for the variational free energy, 
whose derivatives can be handled by means of symbolic 
manipulation programs. Although solving of the corre- 
sponding saddle point equations is cumbersome (even 
numerically), they greatly simplify in the high temper- 
ature case, i.e., when we consider the disordered solution 
TO — tij = Zij = which is isotropic: Xij = x and yij ~ y. 
In that case we obtain (see Supplementary Material): 

log(l + y-2e) = 2 log{e -y)- \Qg{y - x), (21) 



rr (2) 

rp (2) 



4^ 



(13) 



and T^ip^p = 1. For details on the method see, e.g. Ref. 
jssf . In the next subsections we apply the formalism for 
the models © and dJ). 



A. BP approximation for the square lattice case 

Within the representation of the model given by Hamil- 
tonian the orientational order parameter {Q) is basi- 
cally given by the average "magnetization" 



M 



while the coverage is given by 



^7 



(14) 



(15) 



Then, the one site probability densities can be ex- 
pressed as (see Supplementary Material): 



)(5,) = (1 - 0) + i mS, + {^e- l)Sl (16) 

where wc have assumed translational invariance. Defin- 
ing the two-site correlations 



X = y tanh 



(22) 



/3m = -log(2)-21og(l + y-2^^) 
'1 - 



-hSlog 



2 \og{y ~ x). 



(23) 



Working out Eas. (PT|) - (P^ shows that there is only one 
physically meaningful solution. The equilibrium coverage 
9* is then given by the solution of the implicit equation 



where 



F_{9,a) = 2 



-/3m ^ F_(6l*,a), 



\ + y-{e,a) - 29 



y^{e,a){l-a) 
with a = tanh(/3it;/4) and 



(24) 



(25) 



y^[e,a) 



1 r 

2a 



\-a + 2a6- ^ {I - a + 2a0)^ - ia0^ , 

(26) 

while the equilibrium values of the correlations are given 
by y* = y-{0*,a) and x* = ay*. 

To compute the high temperature nematic susceptibil- 
ity wc add to the Hamiltonian ^ a small external field 
B conjugated to to. Then, at temperatures above the 
critical one we can still assume isotropy in the solution 
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(namely, the solution of the saddle point equations that 
converge to the previous one in the limit _B — > 0), so that 
tij — Zij — z and m <^ 1, z ^ 1. This leads to saddle 
point equations which expanded to the lowest order in 
B' = I3B give (see Supplementary Material) 



3— = 4- 



-^B' + 0{m\z^,mz), (27) 



+ 0{m^ ,z'^,mz) 



(28) 



Then, in the linear response regime m ~ xB' and 
z ~ ujB'. In the limit B' — > 0, x is proportional to 
the nematic susceptibility. From the above equations we 
obtain 



r{0* +x*) 

3x* - 0* ' 
)*(.T* +y*) 

3x* - e* 



(29) 
(30) 



The disordered solution becomes unstable whenever 
3a;* = e*. Replacing this condition into Eqs.(pT |) - (P5|) 
we obtain the critical line: 



27 3a - 1 



4 1 



(31) 



and 



= 3 



1 - a 
1 + 3a' 



(32) 



or equivalently: 



4 arctanh ( ^ 



(33) 



In Fig. [5] we compare the Bethe-Peierls critical line 
with those obtained by other techniques. From Eg. ((33)) 
we obtain the following asymptotic behavior when 6^1: 



0.8' 



BP 

Tavares et at 

MC 

RSRG 




square lattice 



0.0 



0.2 



0.4 



1.0 



e 



FIG. 5: Comparison between the square lattice phase dia- 
gram obtained within the Bethe-Peierls (BP) approximation 
and those obtained by other methods: Real Space Renormal- 
ization Group (RSRG) from Ref. Tavares et al. approx- 
imation from Ref. and MC simulations. 



M 

M ^ 



= - [3 ((5(a„ 1)) - 



(35) 



where the broken symmetry direction cr = 1 is taken 
arbitrarily among the different q oriented states ((7 = 3 
in our case). This is a generalization of usual definition 
for the g-state Potts model. In a disordered state we have 
{5[(Ti, 1)) = 9/q,som = 0, while in an ordered state along 
the (T = 1 direction we will have {5{ai, 1)) = so m = 6. 
A conjugated external field to the order parameter ([55)1 
can be considered by adding to the Hamiltonian (HJ a 
term of the form 



B 



(36) 



1 



2 In (0) ' 



(34) 



which agrees qualitatively with the asymptotic behavior 
of Tavares et al. calculation T* ^ -l/3ln{6) 



The coverage is given by 



M 



(37) 



B. BP approximation for the triangular lattice case 

The magnetization (oricntational order parameter) in 
this case is given by 



As in the square lattice case, we will consider here- 
after only isotropic solutions (valid above the transition 
temperature in the B <C 1 limit). We then define the 
correlations 
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X = {6{a,,0)6{<j,,0)), (38) 

y ^ {5ia,,l)S{a,,l)), (39) 

z = {S{a„2)S{aj,2)) ^ {S{a„3)S{a,,3)) , (40) 

t = {S{a,,0)S{aj,l)). (41) 

The one and two sites reduced densities for the spin vari- 
ables (7i = 0, 1, 2, 3 can be expressed as 



p«(a,) = 5]P.(5(a„a) (42) 

a=0 



The physically meaningful solutions of Eqs. (l44l) - ((50)) 
can be obtained in terms of the implicit equation 

0* = G-{e*), (48) 

where 



G_(6l) = - {a + x{e){3 - a)^ 
a 

~y/3[ax{e)[l ~ xie)] + 3x^{e)]'j , (49) 
with a = 2 + and 



(2)/ 



= Pcr,a'S{ai, a) S{aj,a'). 



(43) 



x{e) 



(l_(?)ll/6 



31/6 e/3M/6 615/6 + (1 _ 61)5/6 



(50) 



Applying the normalization and reducibility conditions 
the coefficients and P^.a' can be expressed in terms 
of the parameters m, 0, x, y, z and t (see Supplementary 
Material for details on the calculation). 

First of all we checked the mean-field approximation 
for this model. Inside the CVM formalism, the mean 
field free energy can be obtained by assuming that the 
probability density function is given by p = pl^"* and 
keeping up to the first order term in the cummulant ex- 
pansion of the entropy (33j . With a simple analysis (not 
shown) we found that the mean field approximation pre- 
dicts a first order transition for any value of fi, in com- 
plete disagreement with the numerical simulations. 

Replacing the reduced densities into Eq. (fT2| we obtain 
the BP free energy in terms of the variational parame- 
ters (m, 0, cc, y, z, t) and the corresponding saddle point 
equations (see Supplementary Material). 

At zero field and high enough temperature we have 
a disordered phase, where all ordered states {a ~ 
1,2,3) become equally probable and therefore m = 
{{S{a^, 1)) = 9/q). Also from the definitions (l39|)-(|40l) we 
have that y — z. With some algebra (see Supplementary 
Material) the saddle point equations for the disordered 
solution reduce to 



{i-e-xf{i-ef 



3e'^'', 



9z.T 



[i-e-xf 



(44) 



(45) 



The equilibrium values for the remaining parameters 
are given by x* = x{9*), t* = (1 - 6* - x*)/3 and 



9x* 



-(1 



X ) 



(51) 



Equation p8|) has always at least two solutions for any 
value of /3 and fi: 0* = 1 {x* = 0) and 9* ^ {x* = 1). 
The first one is the meaningful solution in the limit /i — > 
oo. For large but finite values of /i a third solution with 
9* < 1 and 1 — 6** <C 1 emerges, which decreases with 
fi. In the /i <C 1 we obtain the asymptotic behavior (see 
Supplementary Material) 



-/3m 35 



(2 e/3™/3)6 ^ 



-2/3m 39 



(2-}-e^'"/3)ii^ 



1 g/3to/3 

3 2 + e'3"'/3 ■ 



(52) 



(53) 



(54) 



Notice that in the T — !> cxo (/3 — > 0), the correlation 
z* 1/9, as expected. 

At non zero magnetic field B 1 we proceeded as 
in the square lattice case, by expanding the saddle point 
equations and keeping the lowest order in B. This leads 
to the following expression for the nematic susceptibility 
(see Supplementary Material): 



{1-9-x)^ = -x{29 + x-l-3z), (46) 



t= {l-9^x)/3. 



(47) 



X 



1 



'{9z* 



2 126'* 
which diverges when 



5(9^ 



1)' 



(55) 
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1261* - 5(9z* 



+ 1) =0. 



(56) 



For a given value of the chemical potential /^/w, 
Eqs. (j49| and (|58|) must be solved together in order to 
obtain the critical line T* vs. coverage 9. In FiglH] we 
compare a the critical line obtained by numerically solv- 
ing Eqs. (|49| and (|58p with the MC results. In particular, 
in the limit — > oo, when 6^1 and x — !■ 0, we obtain 
from Eqs.dniD-dSll) that 



X 



1 3 + 4e'^"'/3 



(57) 



4 7 _ 4g^tu/3" 
We see that x > for all temperatures T > Tc, with 



1.0 



0.8 



0.6 



J"* 



0.4 



0.2 



0.0 



0.0 




0.2 



0.4 



0.6 







triangular lattice 



0.8 



1.0 



T^; = ^1^- 0.595647. (58) 

and diverges at that temperature, in a clear signature 
of a second order phase transition. Comparing with the 
square lattice result from the previous subsection, we see 
that the critical temperature at ^? = 1 decreases in the 
triangular lattice in a factor w 0.826, which compares 
well with the MC reduction factor « 0.841. 

For finite, but relatively large values of fi/w, Eqs. ip^ 
and (j58p present only one nontrivial solution, that con- 
verges to the value given by Eq.()58p in the /.t — > oo 
limit. However, for fi < fiQ, with ~ —0.65 (which 
corresponds to 6 ~ 0.73), two new non trivial solutions 
emerge, one with 6' <C 1 and the other with 1 — 9 ^ 1. 
As /Li further decreases, the low coverage solution ap- 
proaches the critical one (i.e., that shown in Fig. 
Finally, both solutions collapse at /ic ~ —0.83 (corre- 
sponding to « 0.24) and disappear for n < fXc- Such 
behavior could be indicative of the presence of a first or- 
der transition at low values oi fj., so that the observed 
secondary instabilities in the susceptibility would corre- 
spond to spinodal lines. In that case, a tricritical point 
somewhere along the calculated transition line should be 
expected. Indeed, a similar behavior has been observed 
within the mean field approximation for the square lattice 
[2^ , which disappears in the improved BP approximation 
as we have shown in the previous subsection. However, 
to check whether there is a change in the transition order 
for the triangular lattice case or the observed behavior is 
just spurious, requires a complete minimization analysis 
of the BP free energy in a multidimensional space (taking 
anisotropic ordered solutions into account) which is be- 
yond the scope of the present work. Whatever the case, 
the calculation presented in Fig. [5] should be regarded as 
a high coverage approximation. 

V. REVISITING THE UNIVERSALITY CLASS 

The purpose of this final section is to revisit a num- 
ber of the issues that have emerged during the course of 



FIG. 6: Comparison between the phase diagram obtained 
within the Bethe-Peierls (BP) approximation and MC sim- 
ulations for the triangular lattice. Inset: comparison between 
the BP critical lines for the square and triangular lattices. 



the discussion about the universality class of the SARRs 
model, Refs. [isl - fisl [2l| . For the square lattice case 
[3^ . at intermediate coverages, it was shown that the 
system under study represents an interesting case where 
the use of different statistical ensembles (canonical or 
grand canonical) leads to different and well-established 
universality classes (o = 1 Potts type or q = 2 Potts 
type, respectively) [l3]- In Ref. [ll], Almarza et al. con- 
cluded that the dependence of the universality class of the 
SARRs model on the statistical ensemble, is very likely 
the result of inadequate use of normal scaling to investi- 
gate the critical properties of the constrained (constant 
density) system. However, to date, no completely sat- 
isfactory explanation has been given on the consistency 
of the FSS behavior, in the canonical ensemble, with the 
critical exponents of the ordinary percolation (i.e. 2D 
Potts q = 1 universality class). 

As in III.B.l, we will address here only the triangular 
lattice case. We expect the same universality class for 
chains on honeycomb lattices (with three allowed orien- 
tations), given that the excluded volume term exhibits 
the same symmetry. The critical behavior of the SAARs 
model on a triangular lattice was recently reinvestigated 
by Almarza et al. [2l|. The authors found that the 
isotropic-nematic phase transition occurring in the sys- 
tem is in the 2D Potts g = 3 universality class. This 
conclusion contrasts with that of a previous study in the 
canonical ensemble [l^ which indicates that the tran- 
sition in triangular (and honeycomb) lattices, at inter- 
mediate density , be longs to the q = 1 Potts universality 
class. In Ref. l2l|, Almarza et al. attributed the dis- 
crepancy to the use of the density as the scaling variable 
in Ref . I l9l| . In addition, Almarza et al. have cited a 
paper [35| in which Fisher showed that fixing the density 
in some models corresponds to introducing a constraint 
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FIG. 7: Data collapsing of the Binder cumulant, Ul vs eL^''", 
with the correlation length exponent of the ordinary percola- 
tion (a), and with the renormalized exponent (b). CE means 
canonical ensemble. 
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FIG. 8: Data collapsing of the Binder cumulant, Ul vs eL^^'^ , 
with the correlation length exponent of the 2D Potts q = 3 
universality class, for: grand canonical ensemble (GCE) sim- 
ulations at intermediate coverage (a), and canonical ensemble 
(CE) simulations at full-coverage (b). 



curves corresponding to the Fig. 2(a), where the control 
parameter is the temperature, provides convincing evi- 
dence that the scahng obtained using 1/^=1 = 4/3 is not 
due to the use of the density as the control parameter, 
as claimed by Almarza et al. However, as would be ex- 
pected due to the proximity of the values considered here 
(1^5=1 and v'), a good data collapse with the renormal- 
ized correlation length exponent can also be obtained 
[Fig. 7(b)]. Hence, unlike what happens in the square 
lattice case. Fisher rcnormalization arguments appear to 
be sufficient in the triangular (honeycomb) lattice case. 

Moreover, to check the data presented by Almarza et 
al. [2l[, MC simulations in the grand canonical ensem- 
ble were carried out using an adsorption-desorption al- 
gorithm. It is important to note that the algorithm used 
here is different from that used by Almarza et al. In the 
grand canonical ensemble, the critical behavior was stud- 
ied at the same point of the phase diagram {9c sa 0.5), fix- 
ing the inverse of the reduced temperature to 1/T* = 4.5, 
and varying the chemical potential fi. Very good collapse 
was obtained with v = 5/6 in the scaling plot of Ul [Fig. 
8(a)], thus corroborating the data of Almarza et al. In 
addition, only in the full-lattice limit (9 = 1), canonical 
MC simulations are able to produce results consistent 
with the 2D Potts q = 3 universality class, [see Fig. 
8(b), which shows the collapse of the cumulant curves 
corresponding to the Fig. 2(b)]. 

Finally, we wish to clarify that: (i) We do not hold that 
the universality class of the SARRs model depends on the 
polydispersity of the rods, as was stated by Almarza et 
al. [l^ in reference to our work [T^] . (ii) We agree with 
Almarza et al. that the universality class of the SARRs 
model, in the square lattice, is that of the 2D Ising model, 
whereas that in the triangular and honeycomb lattices, 
is the same as that of the 2D Potts model with q = 
3. (iii) The strong consistency of the results obtained 
in the canonical ensemble with the critical exponents of 
the ordinary percolation (at intermediate coverage, in the 
three lattices considered), warrants an explanation that 
has not yet been given. 



VI. CONCLUSIONS 



that rcnormalizcs the critical exponents. More precisely, 
Almarza et al. have noted that, for the Potts q ~ 3 uni- 
versality class, the renormalized correlation length expo- 
nent I'' is v' = 1^/(1 — a) = 5/4, which is close to the 
value of v for the q = 1 universality class, i'q=i =4/3, 
reported in Ref. |19j . 

In order to test the argument given by Almarza et al., 
a series of MC simulations have been conducted. As in 
Refs. [1^, [13], the distinction between the two univer- 
sality classes is based on the determination of the value 
of which is clearly different for the two universality 
classes under discussion. Then, the scaling behavior can 
be tested by plotting Ul vs eL^^'^ and looking for data 
collapse. As shown in Fig. 7(a), the collapse of the 



In this paper, the main critical properties of self- 
assembled rigid rods on square, triangular and honey- 
comb lattices have been addressed. The results were ob- 
tained by using Monte Carlo simulations in the canonical 
and grand canonical ensembles, finite-size scaling tech- 
niques and theoretical analysis in the framework of the 
Bcthe-Peierls approximation. 

Several conclusions can be drawn from the present 
work. On the one hand, the equilibrium average rod 
length as a function of concentration was calculated by 
MC simulations. In the case of square lattices, computa- 
tional data were compared with theoretical results from 
Tavares et al. A good qualitative agreement was ob- 
served in the range of coverage from to 0.8. However, 
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the disagreement turns out to be significantly large for 
9 > 0.8. In the case of triangular and honeycomb lattices, 
the dependence of the equilibrium average rod length on 
coverage was reported here for the first time. 

The obtained results for J{6) reveal two interesting ob- 
servations: (i) at intermediate coverage {6 ?s 0.5), the 
value of the average rod length coincides with the mini- 
mum value of k {kmm), which allows the formation of a 
nematic phase for a system of monodisperse straight rigid 
/c-mers adsorbed on two-dimensional lattices (square lat- 
tice, kmin = 7 [TTI - [l3t : triangular lattice, kmi7i = 7 
[13, [ill honeycomb lattice, kmin = 11 [3); s-^^d (ii) 
at low coverage, the three curves show the same tendency, 
independently of the lattice geometry (given the range of 
concentrations studied here, this behavior is more ev- 
ident for square and triangular lattices). This finding 
reinforces the idea that the adsorption process behaves 
as a one-dimcnsional problem in the low-coverage (tem- 
perature) regime: particles adsorb forming chains and an 
equilibrium polymerization transition occurs in the sys- 
tem [H. 

On the other hand, and regarding the phase diagram of 
the SARRs, the complete T-9 critical curves correspond- 
ing to triangular and honeycomb lattices have been ob- 
tained by using MC simulation and FSS analysis. In the 
case of triangular lattices, the present study allowed us 
to corroborate previous results obtained from the behav- 
ior of the adsorption isotherms [23| and, in the case of 
honeycomb lattices, the phase diagram has been reported 
here for the first time. 

The simulation phase diagrams were compared with 
analytical data from the BP approximation. BP results 
confirm the whole scenario that emerges from the MC 
simulations, in a clear improvement respect to the mean- 
field approximation, namely: a) a continuous nature of 
the phase transition for any value of (although in the 
triangular lattice case the BP fails in that feature at low 
coverage, the general trend suggests to be an spurious 
effect of the approximation); b) a consistent reduction 
in the critical temperatures for any value 9 when the 
triangular and square lattice models are compared; c) 
an independence on 9 of the critical curves for different 
lattices at very low values of 9 (see inset of Fig|6]) and d) 
a logarithmic decrease with 9 of the critical curve when 
9^0, in agreement with other analytical approach 

In an earlier study [2^ . in which the critical behav- 



ior of SARRs on the square lattice was addressed, it was 
shown that in the full coverage case {9 = 1.0), the Hamil- 
tonian of the SAARs model maps exactly onto the Ising 
model one {q = 2 Potts model) with coupling constant 
Wising = wsARMs/^- In Contrast, the 2D SARRs model 
on the triangular lattice cannot be mapped on g = 3 
Potts model; as can be easily seen from Eq. It is in- 
teresting to note, that through a simple extension of the 
present calculations the critical temperature T* within 
the BP approximation for the isotropic q = 3 Potts model 
results 3 times that of the SARRs on the triangular lat- 
tice. The corresponding comparison between the critical 
temperatures extracted from MC simulations predicts a 
factor 3,3299 « 10/3, in close agreement with the BP 
result. 

Finally, the problem of the universality class of the 
SAARs model was revisited. Since the case correspond- 
ing to square lattices has been widely discussed in Refs. 
Mil, we focused in the case of triangular lattices (the 
same universality class is expected to hold also for hon- 
eycomb lattices with three allowed orientations). Based 
on the calculation of the correlation length exponent v in 
the canonical and grand canonical ensembles, and using 
the Fisher renormalization scheme, we confirmed previ- 
ous results by Almarza et al. |2l|. Namely, the univer- 
sality class of the SARRs model for triangular and hon- 
eycomb lattices is that of the 2D Potts model with q = 3. 
However, the strong consistency of the results obtained 
in the canonical ensemble with the critical exponents of 
the ordinary percolation (at intermediate coverage, in the 
three lattices considered), warrants an explanation that 
has not yet been given. 
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Supplementary information 

Critical behavior of self-assembled rigid rods on two-dimensional lattices: 
Bethe-Peierls approximation and Monte Carlo simulations 

L. G. Lopez, D. H. Linares, A. J. Ramirez-Pastor, D. A. Stariolo and S. A. Cannas 
Here we provide the details of tlie BP approximation presented in the manuscript. 

I. THE BP APPROXIMATION IN THE SQUARE LATTICE CASE 

Following the method outlined in Ref. I] (see chapter 6) the Bethe-Peierls or two-point cluster free 
energy for the square lattice model can be written as: 



F - TvpH + ksT \ (1 - q,) ^ Tr.pf ^ logpf ^ + ^ TT,,,pf] log 



(2) 

<i,j> 



P^,3 



~7 XI + + + W^3){^yn3) + {Xtj + ytj - Zij - W^J){e^.nJ)] - fl^0^ + 

<i,j> i 

+kBT I (1 - 9) E Tr.p!'^ log + E Tr.,,p!? log f^^ \ (1) 

I i <i,J> I 



where: 



Xij 












Zij 




= Tt,,,S,S]p[^^ 






= Ir,,,S^5fpg 



(2) 
(3) 
(4) 
(5) 



We start calculating the one point and two points reduced densities p[^\Si) and p'f'^-{Si, Sj). 
In full generality, the one point density for the model considered can be written as: 



pl^\s,)^A + BS, + CSf. (6) 

Imposing 



1 



mi = {Si)=TriSip\'\Si) 



the values of the coefBcients A, B and C are found to be: 

Analogously for the two point density, consider the expansion: 

pf]iSu Sj) = a + biSi + bjSj + aS^ + cjS^ + dSiSj + e S^S^ + fijSiS] + fjiSjS^ 
Imposing the reducibility conditions 



the following set of equations is obtained: 



3a + 2ci = {1-Oj) 
3a + 2cj = (l-6»i) 
3bi + 2fij = ^rrii 



3bj + 2fji = -rrij 



3 

3ci + 2e = -9i-l 



3cj + 2e = -9,-1 



Actually, these equations are not all linearly independent, since replacing (|T8)) into (|15l) and ([T9)) into p4)) 
we obtain the same equation (Eq. (|2T|) ). Hence, we can choose the following set of linearly independent 
equations 



3a - ^ e = \-e,-e, (20) 
36, + 2/,, = im, (21) 
36, +2/,, - im, (22) 



2 
3 
2 

With the definitions ^ we find d = Xij/A. Also, 



3c, + 2e = I 0,-1 (23) 



3c,+2e = -9,-1 (24) 



j/y = 4a + 4(c, + cj) + 4e (25) 

z,, = 4b,+Af,j (26) 

u;,, = 46, +4/,i (27) 

From Eqs-dll]), dH]), ([Ml) and (E?]) we obtain 



6i = -^{rrii- Zij) (28) 



6, = i(m, - Wy ) (29) 



3 1 

h = 4 ~ 2 '"^ ^^^^ 

/ji = J - ^ TOj (31) 



From Eqs.dini), ^ and ^ we obtain 



= -^-ly^J + l0^ + ^l (32) 



2 " ■' 2 

'^-ly^J+lo,+0, (33) 



3 



a^y^j + l~{e^ + ej) 



(34) 



9 , 3, 



(35) 



Now, from the symmetries of Hamiltonian (Eq.(3) of the manuscript) we can assume — m, 6i = 9, 
but the correlations Xij, yij and Wij and Zij may depend on orientation along the principal axes of the 
lattice. Therefore 



(36) 



and 



1 3 5^ 



4 



-2/,, +1-30 



— z,-i m 



— m 



S,Sf. (37) 



Writing a^^j = x\\ and a:ij = a:^ for pair correlations along and respectively, and the same for yij, 
Wij and Zij, we can write the variational free energy ([I} as: 



K 



[{x\\ + XA_) + {y\\ + y±_) - (z|| - z±_) - {w\\ - w^)] - hO 



S=0,±1 

- E E Pl'^(^i,52)l0gp(')(5i,52)[(e,.fi2) + (e,.fi2)] 



(38) 



Si=0,±l S2=0,±l 



(2) 

where K = I3w. The last trace has contributions from horizontal and vertical links. Calling p\ 2(51, >S'2) 



^(2),para^^^^ 5*2) for thc horizontal case: 

pi'^'P''™(5.,5,) = (2/11 + 1 -20) + i(m-Z||)5, + i(™-^«||)5,+ 



. 3 5, 



{Sj + 5|) + 



4 



^2;il + l-30 



3 1 

— zii m 

4 II 2 



— wii m 

4 II 2 



SjSl (39) 



and a similar expression for the vertical terms p^i2^'^^^{Si,S2), we arrive at a rather long expression 
for the free energy, which can be conveniently handled by a software for symbolic manipulation like 
Mathematica. 



A. High temperature disordered solution 



In the high temperature phase m = w = z = and correlations are isotropic. Deriving the variational 
free energy ([!]) with respect to 9 we find: 



h = -2Log[l+i/||-26']-2Log[l+2/±-26l]+2Log 



-2Log 



2 2 



3 -Log[l -6]+ Log 



Deriving respect to x\\ 
K 

T 



-^Log 
^Log 



a;i| 

X\\ 

2 + -r + 
4 



lM_5. + 2r-l-M + ^ 



2 
(40) 



(41) 



Deriving respect to x±_ : 
K 

T 



-Log 



-Log 



2-1- 



5^ 

2 2 
2 2 



(42) 



Deriving respect to 



- = Log[l + 




- 29] - 


-2Log 


2^11 

2 2_ 


+ ^Log 


+ ^Log 


2H 




13y,| 


- 561 + 2 ^ 


3yii 








4 




2 



2 - 
4 



50 

y 



II , i3yii 



2 -1- 



32/11 , 50 



(43) 



Deriving respect to y±^ : 



^ = Log[l + y^- 29] - 2Log 



--Log 
2 ^ 



132/j 



2 -1- 



-Log 
2 ^ 

32/± 



2-^ 
4 



50 



132/j 



-50 + 2-1- 



3yj 



50 

y 



(44) 



We see that eqa. (PT|) and (H^ are equal and the same happens with P5)) and (HI)) . This confirms 
that the disordered solution is isotropic with respect to correlations and we can write a;|| — xj_ = x and 

y\\ ^y±=y- 

The previous set of equations can be conveniently reduced to: 



/i = - log(2) - 4 log(l + y - 20) + 3 log(l - 0) - 3 log(0) + 4 log(-y + 0) (45) 
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K 
1 



-log 
+ log 



13y 

4 



2-T + ^-56' + 2 
4 4 

2+4 



3y 50 

-1 -A 

2 2 

3v 50 

-1 ^ H 

2 2 



log(y + a;) - log(y - x) 



(46) 



21og[l + t/ 
+ log 



2; 

2+4 



26*] -4 log 
+ H^-50 







+ log 


2-^ 




" 2 




4 


+ 2 ( 






5e'Y 









13y 

4 



561 + 2 -1 



= log(y - x) + log(j/ + x) + 2 log(l + J/ - 261) - 4 log(6' - y) 
where h = These equations can be further simplified to give: 



3y 56* 

y + y 



(47) 



log(l + 2/ - 20) = 2 log(0 - y) - log(j/ - x) 



(48) 



X = V tanh — 



(49) 



h = - log(2) - 2 log(l + y - 20) + 3 log 
Combining Egs. ipSl) and we obtain: 



1 - I 



2 log(j/ - x) 



(50) 



y±(^) 



2a 



1 - a + 2a0 ± y/{l-a + 2a0)2 - 4a02 



(51) 



where a = tanh(ii'/4). It is easy to see that < j/_ < 1 for any value of < < 1 and for any value 
of a. On the other hand, it can be seen that y+ > 1 for any value of < < 1 when a < 0.5. Hence, 
for temperatures ksT/w > 1/(8 arctanh(0.5)) = 0.22756 the only meaningful solution is y_. Replacing 
Eqs. (|5ip into (|50p we obtain two implicit equations for solving as a function of h and namely 



F±{0,a) 



(52) 
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where 



Now, it can be seen that F_{9,a) decreases monotonically with 9, diverging for — > and 
limg^i F^{9,a) = 0, for any value of a. On the other hand, F^{9,a) diverges both in 9 = and 
9 — 1. From the properties y~{0) — 0, — 1, J/+(0) = (1 — a)/a and y+{l) — l/2a, we see that the 

only roots of Eq. (|52]) that satisfy the correct limits 

6' -> 1 y — > 1 for /i -> oo 
^ 2/ -> for h -> -cx) 

is y_. 

B. Near the transition: susceptibilities and critical lines 

Now suppose that we add a small aligning field i?, coupled to m. Now, because of the symmetry 
breaking, one should consider the whole set of parallel and perpendicular correlations, which should be 
different in the two principal directions of the square lattice. Nevertheless, because of the pair approx- 
imation in (H)), parallel and perpendicular correlations appear independently, i.e. the factor involving 
both kinds of correlations do not mix when computing saddle point equations. Then, the saddle point 
equations for each group are exactly the same, implying that parallel and perpendicular quantities them- 
selves are identical. Then, at least in the Bethe-Peierls approximation, there is no symmetry breaking in 
the correlation functions. 

Then, from the full saddle point equations one finds: 

, /m\ f m — z\ , , ^ 

Sarctanh [^—j = 4arctanh I- 1 - B' (54) 

where B' = f3B. If S' <C 1, we can assume m ^ 1 and z ^ 1 and expand 



7 



3 — = 4 B' + 0{m^, z^, mz) 



— — h 0(m , z , mz) 

x + y 9-y 

To this order of approximation, Egs. pS)) . and ((511)) hold. Then, we can assume m = x^' and 
z = LuB' . In the hmit i?' — > 0, x is proportional to the magnetic susceptibility. Replacing in the above 
equations we have: 



9 e-y 



X - w 



x + y e-y 



Solving these equations we obtain 



where x ,y and 9 are solutions of Eas. p5)) - (|5(I| . The disordered solution becomes unstable whenever 
3a; — 9. Replacing these conditions into Egs. pS)) - (15(11) we obtain the critical line: 



. 27 3a - 1 

e = 

4 1 - a 

Notice that, in the limit h oo, we have a = 1/3 or 



(57) 



tanh 



4 / 3 



which is the critical temperature for the Ising model (square lattice) in the Bethe approximation, as 
expected (fcsTc/w — 0.72135). We also have along the critical line 
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II. THE BP APPROXIMATION IN THE TRIANGULAR LATTICE CASE 



We now consider the diluted g' = 3 anisotropic Potts model 



H = -w ^ ^<5(a„a)<5(a,,a)<5(fy,e,)-^^[l-%„0)]--^[3<5(a„l) + (5(a„0)-l] (60) 

<i^j> (7=1 i i 

The one site reduced matrices for these spin variables ai = 0, 1, 2, 3 can be written as 



PrH^0 = E^.<5(a„a) (61) 

(7=0 

From Eq.(37) of the manuscript we have 



(^(<7i,0))=Po = l-^ (62) 
Prom Eq.(35) of the manuscript we have 



{S{ai,l))=Pi = ^{e + 2m) (63) 
Using the symmetry {6{ai, 2)) = {5{ai, 3)) and the normalization condition = 1 we obtain 



Summarizing: 



P2 = Ps = ^(1 - ^'o - Pi) = - m) 



(64) 



pf\^i) = (1 - ^) 5{<J^,Q) + -{e + 2m) 5{a^, 1) + -{9 - m) {5{<j^,2) + (5(a„ 3)) (65) 
We next consider the two-sites reduced matrices 



Pij = ^P<y,a'5(a,,o)5(oj,a') (66) 

where we have that 



Pa^a' = {5{a,,a)5{a,,a')) = ^„ (67) 
From the reducibihty conditions (fT2)) - ((T3l) we have that 



P-,-' = (68) 

cr = 

assuming isotropy (vahd in the disordered state) we define the correlations 



X = {Sia,,0)S{a,,0))^Po,o (69) 

y = (<5(a„l)5(a„l)) = Pi,i (70) 

z = {Sia,,2)Sia,,2)) = (,5(a„ 3) (5(a,, 3)) = F2.2 = ^3,3 (71) 

t = {d{a,,0)d{a,,l)) = Po,i (72) 

and assuming a symmetry under interchange of states a — 2 and a — 3, and using Eqs. (|68p we obtain 

Po,2 = Po,3^^iPo~x~t) = ^il-e-x~ t) (73) 

Pu2 = Pu3 = l{Pi'y-t)^^ i(0 + 2m)-y-t (74) 

P2,3 = ^2 - Po,2 -Pi,2-z = l{e-m) + l{x + y~l) + t-z (75) 
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The variational BP free energy for the triangular lattice can be written as 



F = TrpH + kBT I -5Y,T^iPi'^ iog p^^ + ^ Tv,jp^]logp^] 



<i,j> 
3 



-W. 



N{y + 2z) - p6N - BNm - bNksT ^ log + SNkeT ^ Pa,a' log Pa,a' 

(7=0 (T,cr' 



F/N = -w{y + 2z) - ^iO-Bm 

{ 

+'ikBT\ xlogx + ylogy + 2zlogz + 2tlogt + 2(1 - 6* - a; - i) log 



-5A;bT {{l-e) log(l -e) + -{e + 2m) log 



-{e + 2m) 



+ -^{^ - to) log 



-{e-m) 



-{i-e-x-t) 











) log 






^TO-,-i) 



1/, 2 

-6»+ -m-y 
+2 ( ^(^-m) + i(a; + y-l)+t-z)log 



Deriving respect to m we obtain: 



+ 



2 1 

-(0-m) + -(x + y-l)+t-z 



B _ 5 /6' + 2m 



Deriving respect to 6 we obtain: 



-2 log 



4(61 - m) + 3(a; + y - 1) + 6i - 6z 
6* + 2m - 3(t + y) 



^+log3 = -^log 



{e + 2m){e-mf 



(1-0)3 

+ log [9 + 2m- 3{t + y)]-3 log(l - - a; - t)} 



+ 2 {21og [4(0 - m) + 3(a; + y - 1) + 6t - 6^] + 



Deriving respect to y we obtain: 



BksT 



log6 = logy + log [4(0 - m) + 3(a; + y - 1) + 6i - 6^] 

-2 log [0 + 2m - 3(f + y)] 



Deriving respect to z we obtain: 



w 



SksT 



log6 = log^ - log [4(0 - m) + 3(a; + y - 1) + 6t - 6z] 
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Deriving respect to t we obtain: 



log(2i) - log(l -9-x-t)~log[9 + 2m-3{t + y)] + log [4(61 ^ m) + 3{x + y - 1) + 6t ~ 6z] = (82) 



Deriving respect to x we obtain: 



log a; + log ( - ) - 2 log(l - 6* - a; - t) + log [4(0 - m) + 3{x + y - 1) + 6t ~ 6z] = (83) 



Combining Eqs.lHO]) and ([HT]) we find 



^, fA{9-m) + 3{x + y-l) + 6t-6z\ ^ (z 
I ^ + 2.n^3(^ + ,) j ^ b 

and combining with the rest of equations we arrive to the following set of independent saddle-point 

equations: 



B 5, /0 + 2m\ , (z 

2fc^ = -3^°H^^J" 



(84) 



log 3 



5, 

-3 log 



{Q + 2m)(B -m) 



2 {3 log [6* + 2m - 3(i + y)\ - log(l -9-x~t)-2 log(2i)} 



(85) 



(86) 



log a; + log ( - ) - 21og(l - 6* - a; - t) + log [4(6* - m) + 3(a; + j/ - 1) + 6t - 6z] = 



(87) 



w 



log 6 = log z - log [4(61 - m) + 3(a; + y - 1) + 6i - 6z] 



(88) 



--(t^)--(i) 



(89) 
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A. B = 0: disordered state 



At zero field and high enough temperature we have a disordered phase, where all ordered states (cr = 
1,2,3) become equally probable and therefore m — {{6{ai, 1)) = 0/q). Also from the definitions (|70p 
((73|) we have that y = z and Po.i = Po,2, so 



3t + x = l-0 

which implies that 

4:{e-m) + 3{x + y-l) + 6t-6z = e + 2m-3it + y) = 2e + x-l-3z 



With these conditions we see that Eas. ((84| . (|86|) and (|89|) are automatically satisfied. The remaining 
equations become 



- -51og -— — +61og — (90) 



ksT '=V3(l-6')y' 1-6*- a; 



K7 

log6 = logz-log(26l + a;-l-3z) (91) 



log a; - log ( I ) - 2 log(l - 6 - x) + \og{29 + a: - 1 - 3z) = (92) 



t=l{l-e-x) (93) 



which combined can be rewritten as 



(l-6l-a;)2 
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(95) 



{1-9- xf ^ {2e + x-l- 3z) 



(96) 



t^^{l-0-x) (97) 



From Ea.(|Ml) we can express a; as a function of 9: 



31/6 g/3M/6 615/6 + (1 _ 61)5/6 

From Eq. (|95p we can express z as a function of 9 and x: 



^ ' (98) 



fiw/3 

z= {l-9-xf (99) 

9x ^ ^ ^ ^ 

and replacing the last equation into Eg.dM]) we get a quadratic equation for 9 in terms of x whose 
solutions, combined with Eq. (|98p provide two transcendental equations for 9: 



9^G±{9) 

where 

G± {9) = ^{a + x{'i-a)± Z[ax{l - x) + ix'^]^ (100) 
where a; is given by Eq. (|98l) and 



a = 2 + e^"'/3 

It can be seen that the equation 9 — G+{9) has no solutions, except for — 1, where G+{9) = G-{9). 
The equation 9 — G-{9) has always at least two solutions for any value of /3 and /i: 9 — 1 {x = {)) and 
9 — Q {x — 1). The first one is the meaningful solution in the limit /i — 00. For large but finite values of 
^ a third solution with 9 < 1 and \ — 9 emerges, which decreases with /i. 

14 



Let's consider the /i ^ 1 case. From Ea. (l98)) we have that 



^ 31/6 



and from Eq. (|100p we have 



Combmmg these results we find: 



1-0-7- (101) 



^~20fi 39 

, ..,„/..Mi (102) 



(2 + e^«'/3)ii 



1 e'^'"/^ 
3 2 + e^'"/3 



B. Near the transition: susceptibilities and critical lines 



We now consider the case of a small external field B' ^ 1, where B' = (3B. We can assume 



y = z + e 



(103) 



3t + x + e-l = 6 

where m = 0{B'), e = 0{B') and S = 0{B'). 

Then, expanding Eas. ([5^ - ((5^ and keeping the lowest order in B' , we obtain the following set of linear 
equations for m, e and 5 
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4(2m-2e-(5) = - (104) 



5 e B' , , 

-^™+- = - (105) 



.= ge (106) 



where 



A = 2e + x- l-3z 



C = 1 



and 9, x and z are the solution of Eqs.dMl)-®- Solving Eqs. ([T(M)) -([lM l) . we finally obtain: 



Hence, the susceptibility is 



which diverges when 



B' 9{9z-x + l) 
^ T l20-5i9z-x + l) ('"^^^ 



20(1 -x~9) , , 

S-B' ^^ (108) 



129 - 5{9z - X + 1) 



„, 69 z , ^ 

^-B' —————— (109) 



126* - 5(9z - a; + 1) 



^(9^-^ + 1) (no) 

2 126l-5(9z-a;+l) ^ ' 
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126l-5(9z-x + l) = (111) 
In particular, in the limit /i — oo, when 6^1 and a; — > 0, the susceptibility becomes 



_13 + 4ef«V3 



where we have used Ea. (|103p . 



^ T. Tanaka, Methods m Statistical Physics, Cambridge University Press (2002). 
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